source("../src/communitySim.R")
source("../src/simCancer.R")
source("../src/caseNumbers.R")
if(F) {
# how I got the 0.4899 figure
bob = function(qq) abs(qnorm(0.95, 0, sqrt(qq)) - log(10)/2)
optim(0.5, bob)
}

totalVar = 0.4899106
totalSD = sqrt(totalVar)
propGroupVar = 0.2

doSims=F


parameters = list(
	   geno=log(1.5), env=log(1.25), "geno:env"=log(2),
		 size=150000, probGeno=0.2, 
		 probEnv=0.3, genoMiss= 0.1,
		 envMiss = 0.1, oneEnvPerCommunity =F,
     cancerMissRate=0.15, falseCancerRate=0.00045, 
		 sdGroup= sqrt(propGroupVar) * totalSD, 
		 sdPerson=sqrt(1-propGroupVar)*totalSD, 
		 Enrollmentprobs = c(20000, 40000, 50000, 40000), 
		 agerange = c(35, 69),
		 Followup = c(5, 10, 20, 30),
		 Ncommunity =  50,
		 EqualCommunity=F
)
sig = 0.001
 
populationData = getPopData()  

deathFile = list(
  "F"=read.table("../data/MortalityNoCardiFemale.txt", header=T), 
  "M"=read.table("../data/MortalityNoCardiMale.txt", header=T)
)

DiffCancer= list(
	"F"=lambdaDiffCVD(myfile="../data/IncidenceFemaleCVD.csv",
	 	agegroup),
   M=lambdaDiffCVD(myfile="../data/IncidenceMaleCVD.csv",
	 	agegroup)
	 )	
# common cancer means specific disease being selected. For cvd, we choose all of them  
CommonCancer=NULL

lambda =list(F=list(x=seq(35, 65, by=5), y= exp((-0.5)*0.585)*c(271.95, 343.05,463.46,658.36,798.68,1065.21, 1315.74) / 100000)  ,
M=list(x=seq(35,65, by=5), y=exp((-0.5)*0.585)*c(265.36,407.74,603.52, 846.10, 1077.14, 1407.05, 1649.21)/ 100000) )

propGroupVarseq = c(0.2, 0.35,  0.5, 0.65, 0.8,  0.95)
sdGroup= sqrt(propGroupVarseq) * totalSD 
sdPerson=sqrt(1-propGroupVarseq) * totalSD

SnumberCommunity = c(15, 30 ,50, 80)

Scolour = c("black", "black", "grey", "black")

thecolours =  apply(col2rgb(Scolour)/255, 2, paste, collapse=",")  
thelwd=c(1,1,2,2)
thelty = c(1,2,1,3)

thepch=1 

Nsim = 250


\section{Results}
\label{sec:results}

#\subsection{Incidence Numbers}

#<<caseNumbers,fig=false,results=tex>>=  
if(F){
if(doSims) {
source("../src/caseNumbers.R")

          CommonCancer=NULL

caseNumbers = CaseCI(parameters, lambda, populationData, deathFile, DiffCancer, CommonCancer,verbose=F,
theGender=c("F","M"), Nsim, Percentile=c(0.025, 0.975)) 
  save(caseNumbers, file='caseNumbersCVD.RData')

} else {
	load('caseNumbersCVD.RData')
}
library(abind)
cases = abind(expected=caseNumbers$cases, caseNumbers$CI, along=1)
cases = aperm(cases, c(3,1,4,2))
cases =format(round(cases), width=4,justify='right')
cases=apply(cases, c(1,3,4), function(qq) paste(qq[1], '(', qq[2], ',', qq[3], ')',  collapse=''))
caseMat = matrix(cases, nrow=dim(cases)[[1]],
	dimnames = list(dimnames(cases)[[1]],
		NULL))
		
library(Hmisc)
caseMat = cbind(years=rownames(caseMat), caseMat)
latex(caseMat, file="",
	cgroup =c('Years', 'Colon', 'Stomach'), n.cgroup=c(1,2,2), 
	rowname=NULL, rowlabel=NULL,
#	rowlabel='years',
	colheads=c(' ', rep(dimnames(cases)[[2]],2)),
	caption='Average incidence numbers (with 95\\% prediction intervals) for males and females after various followup periods.',
	label='tab:caseNumbers', where='H', caption.loc='bottom'
)
@

}


\subsection{Relative risk}

<<genoEnvRR,fig=false>>=
if(doSims) {

 

genoEnvRR = seqPowerList(list("geno:env"=log(c(1.5, 1.75, 2, 3, 4))),
	parameters=parameters, CommonCancer=NULL, 
	SimulationTime=Nsim, lambda=lambda,
	populationData=populationData, deathFile) 

save(genoEnvRR, file="genoEnvRRcvd.RData")

} else {
load("genoEnvRRcvd.RData")
}

@
\begin{figure}[H]
\begin{center}
	\subfigure[stroke]{
<<strokeRR>>=

toplot = genoEnvRR[,,'geno:env','stroke',as.character(sig)]



theLabels = rownames(toplot)
theLabels = gsub("[[:alnum:]|:|.]+=", "", theLabels)
theLabels = gsub(",", "", theLabels)
theLabels = signif(exp(as.numeric(theLabels)), 2)

 
  matplot(theLabels, toplot, col=Scolour, xlim=c(1,max(theLabels)),
  lwd=thelwd, lty=thelty,
  type="o",  xlab="relative risk", ylab="power", ylim=c(0,1),
  pch=thepch)

thecolours = toString(paste(dimnames(toplot)[[2]], 
c("( \\\\protect\\\\rule[3pt]{10pt}{1pt} )",
"( - - - )",
"( \\\\textcolor{Gray}{\\\\protect\\\\rule[3pt]{10pt}{2pt}} )",
"( \\\\ldots ) "
) )        )

@	
	}
	\subfigure[hf]{
<<hfRR>>=
toplot = genoEnvRR[,,'geno:env',"hf",as.character(sig)]  
 matplot(theLabels, toplot, col=Scolour, xlim=c(1,max(theLabels)),
  lwd=thelwd, lty=thelty,
  type="o",  xlab="relative risk", ylab="power", ylim=c(0,1),
  pch=thepch)  
@
}
	\subfigure[ihd]{
<<ihdRR>>=
toplot = genoEnvRR[,,'geno:env',"ihd",as.character(sig)]  
 matplot(theLabels, toplot, col=Scolour, xlim=c(1,max(theLabels)),
  lwd=thelwd, lty=thelty,
  type="o",  xlab="relative risk", ylab="power", ylim=c(0,1),
  pch=thepch)  
@
}
	\subfigure[others]{
<<othersRR>>=
toplot = genoEnvRR[,,'geno:env',"others",as.character(sig)]  
 matplot(theLabels, toplot, col=Scolour, xlim=c(1,max(theLabels)),
  lwd=thelwd, lty=thelty,
  type="o",  xlab="relative risk", ylab="power", ylim=c(0,1),
  pch=thepch)  
@
}
	\subfigure[All]{
<<AllRR>>=
toplot = genoEnvRR[,,'geno:env',"All",as.character(sig)]  
 matplot(theLabels, toplot, col=Scolour, xlim=c(1,max(theLabels)),
  lwd=thelwd, lty=thelty,
  type="o",  xlab="relative risk", ylab="power", ylim=c(0,1),
  pch=thepch)  
@
}
	\label{fig:effectSize}
	\end{center}
\end{figure}


\subsection{Prevalence and Misclassification}

<<others,fig=false>>=
SprobGeno =  c(0.05, 0.1, 0.2, 0.4)
SenvMiss = c( 0, 0.05, 0.1, 0.2)

if(doSims) {

probGeno = seqPowerList(list("probGeno"=SprobGeno),
	parameters=parameters, CommonCancer=NULL, 
	SimulationTime=Nsim, lambda=lambda,
	populationData=populationData, deathFile) 

envMiss = seqPowerList(list("envMiss"=SenvMiss ),
	parameters=parameters, CommonCancer=NULL, 
	SimulationTime=Nsim, lambda=lambda,
	populationData=populationData, deathFile) 



save(probGeno, envMiss, file="otherscvd.RData")

} else {
load("otherscvd.RData")
}
@


\begin{figure}[H]
\begin{center}
	\subfigure[Genotype prevalence]{
<<probGeno>>=
Dcancer = 'All'

toplot = probGeno[,,'geno:env',Dcancer,as.character(sig)]
 
matplot(SprobGeno, toplot, col=Scolour,
  lwd=thelwd, lty=thelty,
  type="o",  xlab="proportion", ylab="power", ylim=c(0,1),
  pch=thepch)


@	
	}
	\subfigure[Environmental misclassification rate]{
<<envMiss>>=

toplot = envMiss[,,'geno:env',Dcancer,as.character(sig)]
 
matplot(SenvMiss, toplot, col=Scolour,
  lwd=thelwd, lty=thelty,
  type="o",  xlab="probability", ylab="power", ylim=c(0,1),
  pch=thepch)
 
@
}
	
	\label{fig:others4}
	\end{center}
\end{figure}


\subsection{Community effects and power}

<<numberCommunity,fig=false>>=


parametersNumberCommunity = parameters

parametersNumberCommunity$EqualCommunity=T
parametersNumberCommunity$Followup = 30
parametersNumberCommunity$oneEnvPerCommunity =T


if(doSims) {

         
numbCommBalanced = seqPowerList(list("Ncommunity"=SnumberCommunity),
	parameters=parametersNumberCommunity, CommonCancer=NULL, 
	SimulationTime=Nsim, lambda=lambda,
	populationData=populationData, deathData, DiffCancer) 

parametersNumberCommunity$EqualCommunity=F

 numbCommUnbalanced = seqPowerList(list("Ncommunity"=SnumberCommunity),
	parameters=parameters, CommonCancer=NULL, 
	SimulationTime=Nsim, lambda=lambda,
	populationData=populationData, deathData, DiffCancer) 


save(numbCommBalanced, numbCommUnbalanced, file='numbCommcvd.RData')


parametersNumberCommunity = parameters


parametersNumberCommunity$EqualCommunity=T
parametersNumberCommunity$Followup =  30
parametersNumberCommunity$oneEnvPerCommunity =T

propCommunityVar = seqPowerList(	
	varying =list(sdGroup = sdGroup, sdPerson = sdPerson),
	parameters=parametersNumberCommunity, CommonCancer=NULL, 
	SimulationTime=Nsim, lambda=lambda,
	populationData=populationData, deathFile) 
 
parametersNumberCommunity$EqualCommunity=F

propCommunityVarUnbalanced = seqPowerList(	
	varying =list(sdGroup = sdGroup, sdPerson = sdPerson),
	parameters=parametersNumberCommunity, CommonCancer=NULL, 
	SimulationTime=Nsim, lambda=lambda,
	populationData=populationData, deathFile) 

save(propCommunityVar, propCommunityVarUnbalanced, file="numbCommVarcvd.RData")

} else {
	load("numbCommcvd.RData")
	load("numbCommVarcvd.RData")
}

@



\begin{figure}[H]
\begin{center}
	\subfigure[Changing communities]{
<<commProp>>=
Dfollowup=as.character(parametersNumberCommunity$Followup)
 #$
Dcancer='All'

toplot = cbind(balanced=numbCommBalanced[,Dfollowup,'env',Dcancer,as.character(sig)],
	unbalanced= numbCommUnbalanced[,Dfollowup,'env',Dcancer,as.character(sig)] )

SnumberCommunity = as.integer(gsub("^[[:alnum:]]*=", "", rownames(toplot)))

 
matplot(SnumberCommunity, toplot, col='black',
  lwd=1, lty=1:2,
  type="o",  xlab="number of communities", ylab="power", ylim=c(0,1),
  pch=thepch)
 
@
}
	\subfigure[Changing $\sigma^2/ (\sigma^2+\tau^2)$]{
<<commPropVar>>=

toplot = cbind(balanced=propCommunityVar[,Dfollowup,'env',Dcancer,as.character(sig)],
	unbalanced= propCommunityVarUnbalanced[,Dfollowup,'env',Dcancer,as.character(sig)] )
 
varSeq = matrix(as.numeric(gsub("^[[:alnum:]]+=", "", unlist(strsplit(rownames(toplot), ",")))), ncol=2, byrow=T)
varSeq = (varSeq[,1]^2/(varSeq[,1]^2 + varSeq[,2]^2))

matplot(varSeq, toplot, lty=c(1,2),
  type="o",  xlab="Fraction community variance", ylab="power", ylim=c(0,1), 
  pch=1, col='black')
 
@

}
	\caption{Power for detecting a community-level environmental effect on colorectal cancer after \Sexpr{Dfollowup} years of followup, 
  with a balanced ( --- ) or unbalanced ( - - - ) allocation of individuals to communities, 
  as a function of the number of communities. Significance level is \Sexpr{sig}.}
	\label{fig:communityPower}
	\end{center}
\end{figure}





\subsection{Community effects and bias}



<<biasNcomm,fig=false>>=

paramsNoEffect = parameters
paramsNoEffect$Followup = 30
paramsNoEffect$geno = paramsNoEffect$env = paramsNoEffect[['geno:env']] = 0
paramsNoEffect$oneEnvPerCommunity=T

if(doSims) {
#$
biasCommNumIndiv = seqPowerList(varying =list("Ncommunity"=SnumberCommunity),
	parameters=paramsNoEffect, CommonCancer=NULL,
	SimulationTime=Nsim, lambda=lambda,
	populationData=populationData, deathFile)


 paramsNoEffect$sdGroup =parameters$sdPerson
 paramsNoEffect$sdPerson =parameters$sdGroup

biasCommNumGroup = seqPowerList(list("Ncommunity"=SnumberCommunity),
	parameters=paramsNoEffect, CommonCancer=NULL,
	SimulationTime=Nsim, lambda=lambda,
	populationData=populationData, deathFile)



save(biasCommNumIndiv, biasCommNumGroup,  file="biasNumcvd.RData")

} else {
load("biasNumcvd.RData")
}

@

\begin{figure}[H]
\begin{center}
	\subfigure[20\% community-level variance]{
<<biasCommLow>>=
Dcancer = c('stroke','All')
ltyCancer = c(1,2)
Dfollowup = as.character(paramsNoEffect$Followup )
theylim = c(0, .5)

toplot = biasCommNumIndiv[,Dfollowup,'env',Dcancer,as.character(sig)]
 SnumberCommunity = as.integer(gsub("^[[:alnum:]]*=", "", rownames(toplot)))

matplot(SnumberCommunity, toplot,
  lty=ltyCancer,
  type="o",  xlab="communities", ylab="Type 1 error", ylim=theylim,
  pch=thepch, col='black')

abline(h=sig, lty=3)
@
	}
	\subfigure[80\% community-level variance]{
<<biasCommHigh>>=

toplot = biasCommNumGroup[,Dfollowup,'env',Dcancer,as.character(sig)]

matplot(SnumberCommunity, toplot,
  lty=ltyCancer,
  type="o",  xlab="communities", ylab="Type 1 error",   ylim=theylim,
  pch=thepch, col='black')

abline(h=sig, lty=3)
@
}
	\caption{Proportion of type 1 errors for detecting a community-level environment effect, for \Sexpr{dimnames(biasCommNumGroup)[[4]][1]} ( --- ) and \Sexpr{dimnames(biasCommNumGroup)[[4]][2]} ( - - - ) cancer at \Sexpr{Dfollowup} years followup.
    Significance level is \Sexpr{sig} \dottedline.}
	\label{fig:biasNumComm}
	\end{center}
\end{figure}



<<biasPropSDgroup,fig=false>>=
paramsNoEffect = parameters
paramsNoEffect$Followup = 30
paramsNoEffect$geno = paramsNoEffect$env = paramsNoEffect[['geno:env']] = 0
paramsNoEffect$oneEnvPerCommunity=T

sig2='0.01'

if(doSims) {


biasComm = seqPowerList(varying =list(sdGroup = sdGroup, sdPerson = sdPerson),
	parameters=paramsNoEffect, CommonCancer=NULL, 
	SimulationTime=Nsim, lambda=lambda,
	populationData=populationData, deathData, DiffCancer) 

paramsNoEffect$oneEnvPerCommunity=F
#$
biasIndiv = seqPowerList(	varying =list(sdGroup = sdGroup, sdPerson = sdPerson),
	parameters=paramsNoEffect, CommonCancer=NULL, 
	SimulationTime=Nsim, lambda=lambda,
	populationData=populationData, deathData, DiffCancer) 
 


save(biasComm, biasIndiv,  file="biascvd.RData")
} else {
load("biascvd.RData")
}
@



\begin{figure}[H]
\begin{center}
	\subfigure[Community-level effect, \Sexpr{sig} significance]{
<<biasComm>>=
Dcancer = c('stroke','All')
Dfollowup = as.character(paramsNoEffect$Followup )
 #$
toplot = biasComm[,Dfollowup,'env',Dcancer,as.character(sig)]
varSeq = matrix(as.numeric(gsub("^[[:alnum:]]+=", "", unlist(strsplit(rownames(toplot), ",")))), ncol=2, byrow=T)
varSeq = (varSeq[,1]^2/(varSeq[,1]^2 + varSeq[,2]^2))

matplot(varSeq, toplot, col=Scolour,
  lwd=thelwd, lty=thelty,
  type="o",  xlab="Fraction community variance", ylab="Type I error",
  pch=thepch, ylim=c(0, 0.02))

abline(h=sig, lty=3)
@	
	}
	\subfigure[Individual-level effect, \Sexpr{sig2} significance]{
<<biasIndiv>>=

toplot = biasIndiv[,,'env',Dcancer,as.character(sig2)]
 
matplot(varSeq, toplot, col=Scolour,
  lwd=thelwd, lty=thelty,
  type="o",  xlab="Fraction community variance", ylab="Type I error", 
  pch=thepch, ylim=c(0, 0.02))
 abline(h=sig2, lty=3)
@
}
	\caption{Proportion of type 1 errors  for detecting an environment effect, for \Sexpr{dimnames(biasComm)[[4]][1]} ( --- ) and \Sexpr{dimnames(biasComm)[[4]][2]} ( - - - ) cancer, with significance level \dottedline. Follow up times of \Sexpr{Dfollowup}  years.  }
	\label{fig:biasFracVar}
	\end{center}
\end{figure}

